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Abstract. Using atomistic computer simulations, we study how ion irradiation can 
be used to alter the morphology of a graphene monolayer by introducing defects of 
specific type, and to cut graphene sheets. Based on the results of our analytical 
potential molecular dynamics simulations, a kinetic Monte Carlo code is developed 
for modelling morphological changes in a graphene monolayer under irradiation at 
macroscopic time scales. Impacts of He, Ne, Ar, Kr, Xe and Ga ions with kinetic 
energies ranging from tens of eV to 10 MeV and angles of incidence between 0° and 
88° are studied. Our results provide microscopic insights into the response of graphene 
to ion irradiation and can directly be used for the optimization of graphene cutting 
and patterning with focused ion beams. 
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1. Introduction 

Recent reports on large-scale production of graphene [H El E], including growth 
of centimeter-size sheets on copper surfaces [4j, have brought closer the utilization 
of graphene's excellent electronic properties [5] in future electronic devices. In 
particular, quantum dots |6], nanoribbons [7], and antidot lattices |8l|9], which provide 
electron confinement within the graphene plane, have received considerable attention. 
Production of such structures is based on graphene patterning by various techniques such 
as electron beam lithography combined with reactive ion etching [6l |TOl [11], chemical 
methods including unfolding of carbon nanotubes fl2\ [13] and graphene cutting with a 
focused electron beam [14j . 

Since focused ion beams (FIB) are already routinely employed in today's 
semiconductor industry, this method could become an alternative approach. Indeed, 
cutting and patterning graphene using a FIB with a high spatial resolution was recently 
demonstrated [IHl UHl [17]. The method is based on using a 30 keV He ion beam in a 
helium ion microscope to directly sputter carbon atoms from predetermined areas of 
graphene sheets, either suspended or deposited on a substrate. Etched gaps down to 
10 nm have been reported, with sharp edges and little evident damage or doping to the 
sample. Concurrently with the experiments on He beam-based patterning of graphene, 
the possibility of using Ar ion irradiation has also been studied [H]. 

Efficient use of ion beams and optimization of the graphene cutting process require 
detailed microscopic knowledge of damage production mechanisms and types of defects 
created by the energetic ions in the sample. In order to get insight into the cutting 
process, the interaction of He ions with the target was modeled in Ref. [16] by a 
semi-empirical method [19] based on the binary-collision approximation, combined with 
statistical algorithms to calculate how a moving ion transfers its energy to the target 
atoms. This approach implemented in the TRIM software package [20] gives reasonable 
results for bulk materials. However, as has been pointed out [211 l22], the theory of 
irradiation effects developed for bulk materials does not always work at the nanoscale. 
In particular, it was recently demonstrated [23] that this approach cannot be applied to 
graphene, as the sample is treated as an amorphous matrix with a homogeneous mass 
density while the explicit account for the atomic structure is required for atomically 
thin targets. 

In this study, we use analytical potential (AP) molecular dynamics (MD) 
simulations, a much more accurate method than the one implemented in the TRIM 
code, to study defect production in graphene under ion irradiation. Our ultimate goal 
is to provide the means for determining optimum parameters, such as ion mass, energy 
and angle of incidence for graphene cutting, which would enable the production of 
smooth edges with the minimum number of defects at faster cutting rates. We simulated 
more than two million impacts of energetic ions onto suspended graphene and gathered 
statistics on types and abundances of defects for He, Ne, Ar, Kr, Xe and Ga ions with 
energies ranging from tens of eV to 10 MeV. The role of the angle of incidence of the 
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ions was studied in detail. To establish a direct link to the experiments, we further 
developed a kinetic Monte Carlo (kMC) code [2l], which utilizes the statistics from 
the MD simulations for predicting the evolution of graphene under ion irradiation at 
macroscopic timescales. This allows modelling the behavior of irradiated graphene under 
realistic experimental conditions in, e.g., FIB systems. 

2. Methods 

In our MD simulations (a general introduction to this method can be found in Ref. [25] ) 
the carbon-carbon interaction in graphene was modelled using the second-generation 
reactive empirical bond-order Brenner potential [26]. The bond conjugation term, which 
is not expected to be important in irradiation processes [27] , was omitted. This approach 
has previously been successfully used in modelling the irradiation response of graphene 
and other carbon nanostructures [23l EHl EHl [301 El] . The interaction between energetic 
noble gas ions (He, Ne, Ar, Kr and Xe) and target carbon atoms was simulated using the 
Ziegler-Biersack-Littmark universal repulsive potential [19] . In addition to the noble gas 
ions, irradiation with Ga ions was simulated using the same universal repulsive potential, 
as gallium is a typical ion species used in FIB systems. Chemical effects are expected to 
be important only at low ion energies (< 250 eV) and they should not play any role at 
typical operating energies (~30 keV) of FIB systems, especially in such a thin target, so 
that the use of the purely repulsive potential is justified also for the Ga-C interaction. 

For a direct analogy with the experiments, we use the term "ion" throughout this 
article, although the charge of the incoming ion is not explicitly considered, as this is 
beyond the AP-MD method, and more importantly, the effects of low charge states are 
negligible [32]. The AP-MD approach is computationally efficient enough for running 
the massive number of simulations (more than two million runs in total) required to get 
statistically meaningful results in the large parameter space explored. 

It is well known that electronic stopping is the main mechanism of energy 
transfer from an ion to any solid target at high ion energies [19]. At the same time, 
experiments [331 El] indicate that energy deposited into electronic degrees of freedom 
of graphite gives rise to defects only if electronic stopping power exceeds a critical 
value of 7 keV/nm. Taking into account the excellent charge and heat conductance of 
graphene, similar behaviour can be expected for this material [32]. As the typical 
electronic stopping value for all the ion/energy combinations used in our study is 
less than 0.7 keV/nm (we calculated electronic stopping power using the approach of 
Ziegler, Biersack, and Littmark [19]), and even the highest values for high-energy Xe 
(~ 4.5 keV/nm) are well below the critical one, electronic stopping effects were not 
taken into account. 

The simulated graphene target consisted of 800 carbon atoms. Since thermal 
excitations due to temperature-related atomic motion do not play a significant role 
in the momentum transfer between the impinging ion and target atoms (typical phonon 
energies are within tens of meV per atom), the initial target temperature was chosen 
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to be K. After ion impacts the system was cooled down using the Berendsen 
thermostat [35] at the edges, with At = 10 fs as the time constant. Such a setup was 
used to model the dissipation of heat generated by the ion impact as would happen if the 
target was part of a larger graphene sheet. Additional simulations with the thermostat 
turned off showed, however, that the results are not sensitive to its parametrization, as 
the damage creation processes typically take place before the heat wave generated by 
the ion impact even reaches the system edges, barring a few cases in the high ion energy 
and high angle of incidence regime. Adaptive time step MD code parcas [36] was used 
for the simulations. Within the code, the simulation time step is dynamically adjusted 
based on the fastest moving particles in the system. This resulted in time steps ranging 
from the attosecond scale to close to one femtosecond. 

The system was allowed to relax for one picosecond after the ion impact. During 
this time the system had typically reached a local energy minimum. However, any 
reasonable simulation time is too short for the system to always find the most stable 
local configuration. To address this, the resulting structures were annealed for another 
picosecond at 1500 K and eventually cooled down to K to facilitate further relaxation 
of the system and removal of at least the most unstable atomic configurations. Over 
macroscopic time scales at experimentally relevant temperatures (room temperature and 
above) the created defects are expected to reconstruct into patches of non-hexagonal 
carbon rings similar to those presented in Ref. [37] . Although the AP-MD method 
cannot capture the reconstruction processes, the sizes of these patches are defined by 
the initial damage which is statistically predicted by our results. 

Impact points for the ions were chosen randomly within the minimum irreducible 
area of the graphene lattice. The direction of the incoming ion was determined by 
randomly selecting an inclination angle 9 G [0,88°] and azimuthal angle G [0,360°[ 
independently so that for each 6 an average result could be generated over tp. A 
schematic illustration of the simulation setup is presented in Fig.[T]along with an example 
of the effect of changing the 9 angle. 

A link between the theoretical data on damage production and experiments can be 
provided within the framework of the kMC method. We have previously implemented 
such a model for comparing theoretical results to experiments in the cases of electron 
irradiation of carbon nanotubes [38] and hexagonal boron nitride monolayers [39]. Here, 
we developed a kMC model describing the response of graphene to ion irradiation by 
using the results from our AP-MD simulations as the source of damage event rates. 
Within this model, the impact event rate is derived from ion beam current, and the 
rates of producing any type of defects are derived by multiplying the ion impact rate 
by defect creation probability for the selected combination of ion species, irradiation 
energy and angle of incidence. We employed data on the defect size (area initially 
covered by perfect hexagonal carbon rings transformed into other polygons, termed 
"amorphized area" within our kMC model) and the number of sputtered carbon atoms 
for each parameter combination. Because the MD simulations were always carried out 
on a pristine target, one should be careful when interpreting the kMC results when the 




Figure 1. Ball-and-stick presentation of the simulation setup, (a) Ion impact in 
direction perpendicular to the graphene plane, (b) Structure of the graphene target 
after the impact (complex vacancy) , where the damage is caused by an in-plane collision 
cascade, (c) Ion impact at an angle 9 = 30° toward exactly the same spot as in (a), 
(d) Structure after the impact (single vacancy), where the C atom is recoiled out of 
the plane, thus creating no further damage. The crosses mark the ion impact points 
and the squares the original locations of sputtered atoms. Both impacts were carried 
out with 2 MeV Ar+ ions. 

defects start to overlap. However, as it is easier to displace under-coordinated atoms, 
our kMC results can be used to estimate the lower limit for the number of sputtered 
atoms even for overlapping defects. 

3. Results and Discussion 

3.1. Dynamical simulations of defect production 

When analyzing the results of the MD simulations, we categorized the defect structures 
produced in graphene into single vacancies (one missing atom from otherwise intact 
structure), double vacancies (two missing atoms and no further damage), complex 
vacancies (defect structures with missing atoms other than single and double vacancies) 
and amorphous regions (defect structures with no missing atoms, the simplest example 
of such a structure being the Stone- Wales defect [10]). Experimental examples of these 
structures are presented in a recent review article [H], except for complex vacancies, 
examples of which are presented in Ref. [37]. We calculated the probabilities of 
producing each of these defect types for all ion species and the whole range of energies 
(35 eV - 10 MeV) and angles of incidence {6 G [0,88°]) considered, as presented in 
Fig. [21 Also the probabilities of creating any defect was calculated. 

As is evident from Fig. |2l varying any of the parameters has a drastic effect on 
the produced defect types and their abundances. At energies in the keV range, single 
vacancies are the typical defects when the angle of incidence is perpendicular to the 
graphene sheet. Upon tilting the ion beam, first double vacancies become the dominant 
defect type after which complex vacancies become the most common type of defect. This 
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can be attributed to the increased projected density of graphene, when viewed from a 
grazing angle. The maxima in single vacancy production are shifted towards higher ion 
energies as the angle of incidence is increased, which can be useful when ion irradiation 
with a fixed energy is used for cutting graphene. The location of the maxima in complex 
vacancy production first moves toward lower energies with increasing angle of incidence 
after which it again moves to higher energies. Therefore, if graphene is to be used as a 
membrane in external beam experiments as suggested in Ref. |23], the resilience of the 
membrane can be improved by slightly tilting it. For example, in the case of 2 MeV 
Ne ions, tilting the graphene sheet by ~ 18° will result in a twofold decrease in the 
sputtering rate, as will be discussed in more detail later in this article. 

A range of parameters where no defects are created corresponds to the large angle, 
low energy region of the graphs in Fig. |2l Low defect production has the same origin as 
under ion channeling conditions in crystalline solids. When coming in at a grazing angle, 
the ion interacts with a long row of atoms, and if the energy is low enough, depending 
on the angle of incidence, none of the target atoms receives enough momentum to 
be displaced. Instead, the ion is reflected away from the graphene sheet. However, 
if ion energy is increased, the ion will penetrate the sheet, and it typically creates 
significant damage along the way, as can be seen from the steeply rising complex vacancy 
probabilities in the low energy-large angle regimes in Fig. |2l Amorphization events are 
observed mostly at low ion energies (<~ 1 keV), where the ions have barely enough 
kinetic energy to displace a target atom, but the displaced atom remains bonded to the 
sheet as an adatom or a part of a local amorphization. 

3.2. Statistical model for continuous defect production 

General trends and dominant defect types for specific irradiation conditions can be 
directly obtained from Fig. O However, making quantitative predictions on the evolution 
of a graphene target under ion irradiation is not straightforward, as both variations in 
the defect sizes and the number of sputtered atoms must be taken into account. Also, 
the probability distribution for creating each defect is not Gaussian, which means that 
mean size and standard deviation are not sufficient for describing the damage caused to 
graphene. For this reason, we developed a kMC code, which directly uses the statistics 
on the sizes of the defects and number of sputtered atoms, extracted from the MD 
simulations. This code can be used for predicting the evolution of a graphene sheet under 
ion irradiation at macroscopic time scales. The code has a web-based interface [21] and 
is available for public use. Both a visual representation of the sample and quantitative 
estimates of the total amorphized area and the number of sputtered atoms can be 
obtained using the program. As the original data is collected for individual ion impacts 
on pristine graphene, we stress that the code underestimates the produced damage at 
high ion doses due to a drop in the displacement threshold energy for under-coordinated 
atoms [22] . 

To demonstrate the potential of the kMC code, in Figs. [3] and H] we present example 
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Figure 2. Probabilities of producing single vacancies, double vacancies, complex 
vacancies, local amorphizations and any defect (any modification to the pristine 
structure) in graphene under ion irradiation as functions of angle of incidence 9 and 
ion energy K. Data for Ga is presented only for energies K >0.5 keV. Data for lower 
energies is excluded as chemical effects were not accounted for in the Ga-C interactions. 
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Ne\ 0.65 keV. 6; 0°, *: 1x10" om"^ Kr*, 0.4 keV, 9: 0°, IxlO" cm-= 
L: 1.3 %, A: 7.7 % L: 3.6 %, A: 16 % 



Ne*, 2 MeV, 9:0°, 1x10'= cm"^ Ne*, 2 MeV. 9: 17.5°, 1);10«cm-= 
L: 0.22 %, A: 3.6 % L: 0.16 %, A: l.g % 



He*. 30 keV. 9: 0°, dose: 0.001 pC He*. 30 keV. 6: 62.5°, dose: 001 pC 
L: 1.9 %, A: 17 % L: 2.4 %, ^: 12.5 % 



Figure 3. Example results of kinetic Monte Carlo simulations of graphene evolution 
under ion irradiation. Blue circles mark single vacancies, whereas the red ones stand for 
other vacancy- type defects, and the green ones for amorphized areas with no missing 
atoms. Percentages of lost carbon atoms (L) and amorphized area (A) are given 
for each case. Note that vacancy-type defects also contribute to the percentage of 
amorphized area. All the numerical values were averaged over 25 simulations. The 
length of the scale bars is 5 nm in all the panels. The pink areas in panels (e) and (f) 
indicate the spot area of a focused ion beam limited at the full width at half maximum 
of the spot intensity. $ is irradiation fiuence as in ion impacts per unit area and dose 
stands for total accumulated charge of the impinging ions. 



cases simulated with the code. If graphene is uniformly irradiated with 650 eV Ne ions at 
normal direction to the sheet, nearly exclusive production of single vacancies is observed 
(Fig. |3^). Changing the ion to 400 eV Kr will result in uniform distribution of single 
and double vacancies with roughly a 1/1 ratio (Fig. EJ^b)). For high energy ions with 
incident direction normal to the sheet the probability of creating any defect is greatly 
reduced, as is evident from Fig. [3t^c) for the case of 2 MeV Ne: Although the fiuence is 
increased by one order of magnitude as compared to the previous examples, the total 
number of defects is lower. The individual defects are much larger, however, which 
can be attributed to the fact that during the interaction of a target atom and the 
fast ion the target atom is practically immobile, resulting in symmetrical momentum 
transfer in the direction of the ion trajectory. This leads to recoiled atoms moving very 
close to a direction perpendicular to the ion trajectory. In the case of ion trajectory 
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Figure 4. Results of kMC simulations of graphene evolution under ion irradiation 
with paranictrization similar to what was used in experiments in Ref. |15| (panel (a)) 
and with optimized irradiation angle with regard to cutting efficiency (panels (b) and 
(c)). Blue circles mark single vacancies, the red ones stand for other vacancy-type 
defects, and the green ones for amorphized areas with no missing atoms. The pink 
area indicates the area swept by the ion beam limited at the full width at half maximum 
of the beam spot. Percentages of lost carbon atoms (L) and amorphized area (A) are 
given for each case. The scale bars are 5 nm. It is evident that by tilting the ion beam, 
the sputtering rate is increased significantly and the creation of local amorphizations 
is avoided. Line dose stands for total accumulated charge through charged ion impacts 
per unit length. 

perpendicular to the graphene sheet, this leads to in-plane collision cascades and large 
defect structures. This effect is also demonstrated in Fig. [H If the direction of the ions 
is tilted away from the normal of the graphene plane, the typical direction of the recoiled 
atoms is correspondingly tilted out of the graphene plane, which results in smaller defect 
structures due to the absence of collision cascades, as can be seen in Fig. |3t^d). Thus 
the resilience of graphene under high energy ion irradiation can be improved by tilting 
it, as was suggested above. 

3.3. Optimizing ion beam cutting of graphene 

As experimentally demonstrated by Lemme et al. [IS] and Pickard et al. [T7] , a focused 
He ion beam can be used to cut graphene with a very high precision. The used 
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acceleration voltage was 30 kV and the cutting was performed with the ion beam 
pointed in the normal direction of the sample surface. The effects of such irradiation on 
graphene are shown in Fig. |3](e) assuming a focused spot with a Gaussian distribution 
with a full width at half maximum (FWHM) of 6 nm. Approximately half of the damage 
events lead to sputtering of target atoms and half of the events lead to amorphizations. 
If clean cutting is to be achieved, the amorphization events are not desirable. The 
situation can be improved by tilting the beam by approximately 60° (Fig. Et^f)). As the 
projected atomic density of the target is now increased, the defect creation probability 
is correspondingly higher. The percentage of sputtered atoms inside the spot area 
(illustrated with a white line along the FWHM edge) is increased with the same total 
dose even though the spot area is larger. Also, as demonstrated previously, tilting the 
beam at high ion energies leads to decreased probability of creating in-plane collision 
cascades and local amorphizations, in favor of single vacancies. This leads ultimately to 
cleaner edges in cuts made with a FIB. 

To further illustrate what can be achieved with parameter optimization when 
cutting graphene with a FIB, three examples of the use of the kMC code in linear 
scan mode are presented in Fig. HI In these examples, the beam spot is moved in the 
vertical direction. The spot width is 10 nm (FWHM) and the ion beam direction is tilted 
towards the direction of the spot movement. The ion energy and species are similar to 
what was used in Ref. [I5]. Line dose of 8 nC/cm was reported to not be adequate for 
making a cut in graphene. This observation is supported by our simulations (Fig. Hl^a)): 
Only ~6 % of the carbon atoms within the cut area are sputtered, although much of 
the area is amorphized. However, when the ion beam is tilted by 62.5° there is almost 
a threefold increase in the number of sputtered atoms (Fig. IHb)). To compare the 
effectiveness of graphene cutting at different angles of incidence, a third simulation was 
conducted, where the line dose was decreased so that the number of sputtered atoms is 
approximately the same as in the first case. Comparing the results in Figs. Hl^a) andlH^c) 
shows that tilting the He beam by ~60° gives a threefold increase in the cutting efficiency 
of the beam, while the amorphized area outside the cut is decreased considerably. 

4. Conclusions 

To conclude, we studied the effects of ion irradiation on graphene using atomistic 
computer simulations. The role of the angle of ion incidence in a wide range of energies 
was investigated, and the types and concentrations of defects were identified for various 
ions. The dramatic effect of the angle of incidence on defect production demonstrates 
the fundamental difference of strictly two-dimensional graphene from traditional bulk 
targets. The peculiarities of graphene's response to ion irradiation can be used to 
gain detailed control over produced defect types and their abundances. The presented 
publicly available computer code [21] enables one to make quantitative predictions of 
defect production in graphene under energetic ion bombardment. This information is 
needed in order to controUably create specific types of defects, or when graphene is to 



Cutting and controlled modification of graphene with ion beams 



11 



be nanomachined by a FIB. To illustrate the latter, we showed an example where ~ 60° 
tilting of the sample gave a threefold increase in sputtering and reduced amorphized 
areas in the sample. Further on, with respect to the possible use of graphene membranes 
in ion-beam analysis [23], we showed that the resilience of the membrane can be improved 
for high energy ions by choosing an optimum angle with regard to the direction of 
the ion beam. Although we looked specifically at noble gas irradiation, our results 
can also provide insights into the response of graphene to irradiation by other species, 
e.g., nitrogen, as in very recent experiments on irradiation mediated nitrogen doping of 
graphene [S]. 
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